home *** CD-ROM | disk | FTP | other *** search
/ IRIX Base Documentation 2002 November / SGI IRIX Base Documentation 2002 November.iso / usr / share / catman / p_man / cat3 / SCSL / dlatps.z / dlatps
Encoding:
Text File  |  2002-10-03  |  8.6 KB  |  265 lines
  1.  
  2.  
  3.  
  4. DDDDLLLLAAAATTTTPPPPSSSS((((3333SSSS))))                                                          DDDDLLLLAAAATTTTPPPPSSSS((((3333SSSS))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      DLATPS - solve one of the triangular systems  A *x = s*b or A'*x = s*b
  10.      with scaling to prevent overflow, where A is an upper or lower triangular
  11.      matrix stored in packed form
  12.  
  13. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  14.      SUBROUTINE DLATPS( UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE, CNORM,
  15.                         INFO )
  16.  
  17.          CHARACTER      DIAG, NORMIN, TRANS, UPLO
  18.  
  19.          INTEGER        INFO, N
  20.  
  21.          DOUBLE         PRECISION SCALE
  22.  
  23.          DOUBLE         PRECISION AP( * ), CNORM( * ), X( * )
  24.  
  25. IIIIMMMMPPPPLLLLEEEEMMMMEEEENNNNTTTTAAAATTTTIIIIOOOONNNN
  26.      These routines are part of the SCSL Scientific Library and can be loaded
  27.      using either the -lscs or the -lscs_mp option.  The -lscs_mp option
  28.      directs the linker to use the multi-processor version of the library.
  29.  
  30.      When linking to SCSL with -lscs or -lscs_mp, the default integer size is
  31.      4 bytes (32 bits). Another version of SCSL is available in which integers
  32.      are 8 bytes (64 bits).  This version allows the user access to larger
  33.      memory sizes and helps when porting legacy Cray codes.  It can be loaded
  34.      by using the -lscs_i8 option or the -lscs_i8_mp option. A program may use
  35.      only one of the two versions; 4-byte integer and 8-byte integer library
  36.      calls cannot be mixed.
  37.  
  38. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  39.      DLATPS solves one of the triangular systems A *x = s*b or A'*x = s*b with
  40.      scaling to prevent overflow, where A is an upper or lower triangular
  41.      matrix stored in packed form. Here A' denotes the transpose of A, x and b
  42.      are n-element vectors, and s is a scaling factor, usually less than or
  43.      equal to 1, chosen so that the components of x will be less than the
  44.      overflow threshold.  If the unscaled problem will not cause overflow, the
  45.      Level 2 BLAS routine DTPSV is called. If the matrix A is singular (A(j,j)
  46.      = 0 for some j), then s is set to 0 and a non-trivial solution to A*x = 0
  47.      is returned.
  48.  
  49.  
  50. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  51.      UPLO    (input) CHARACTER*1
  52.              Specifies whether the matrix A is upper or lower triangular.  =
  53.              'U':  Upper triangular
  54.              = 'L':  Lower triangular
  55.  
  56.      TRANS   (input) CHARACTER*1
  57.              Specifies the operation applied to A.  = 'N':  Solve A * x = s*b
  58.              (No transpose)
  59.              = 'T':  Solve A'* x = s*b  (Transpose)
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. DDDDLLLLAAAATTTTPPPPSSSS((((3333SSSS))))                                                          DDDDLLLLAAAATTTTPPPPSSSS((((3333SSSS))))
  71.  
  72.  
  73.  
  74.              = 'C':  Solve A'* x = s*b  (Conjugate transpose = Transpose)
  75.  
  76.      DIAG    (input) CHARACTER*1
  77.              Specifies whether or not the matrix A is unit triangular.  = 'N':
  78.              Non-unit triangular
  79.              = 'U':  Unit triangular
  80.  
  81.      NORMIN  (input) CHARACTER*1
  82.              Specifies whether CNORM has been set or not.  = 'Y':  CNORM
  83.              contains the column norms on entry
  84.              = 'N':  CNORM is not set on entry.  On exit, the norms will be
  85.              computed and stored in CNORM.
  86.  
  87.      N       (input) INTEGER
  88.              The order of the matrix A.  N >= 0.
  89.  
  90.      AP      (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
  91.              The upper or lower triangular matrix A, packed columnwise in a
  92.              linear array.  The j-th column of A is stored in the array AP as
  93.              follows:  if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  94.              if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
  95.  
  96.      X       (input/output) DOUBLE PRECISION array, dimension (N)
  97.              On entry, the right hand side b of the triangular system.  On
  98.              exit, X is overwritten by the solution vector x.
  99.  
  100.      SCALE   (output) DOUBLE PRECISION
  101.              The scaling factor s for the triangular system A * x = s*b  or
  102.              A'* x = s*b.  If SCALE = 0, the matrix A is singular or badly
  103.              scaled, and the vector x is an exact or approximate solution to
  104.              A*x = 0.
  105.  
  106.      CNORM   (input or output) DOUBLE PRECISION array, dimension (N)
  107.  
  108.              If NORMIN = 'Y', CNORM is an input argument and CNORM(j) contains
  109.              the norm of the off-diagonal part of the j-th column of A.  If
  110.              TRANS = 'N', CNORM(j) must be greater than or equal to the
  111.              infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) must be
  112.              greater than or equal to the 1-norm.
  113.  
  114.              If NORMIN = 'N', CNORM is an output argument and CNORM(j) returns
  115.              the 1-norm of the offdiagonal part of the j-th column of A.
  116.  
  117.      INFO    (output) INTEGER
  118.              = 0:  successful exit
  119.              < 0:  if INFO = -k, the k-th argument had an illegal value
  120.  
  121. FFFFUUUURRRRTTTTHHHHEEEERRRR DDDDEEEETTTTAAAAIIIILLLLSSSS
  122.      A rough bound on x is computed; if that is less than overflow, DTPSV is
  123.      called, otherwise, specific code is used which checks for possible
  124.      overflow or divide-by-zero at every operation.
  125.  
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. DDDDLLLLAAAATTTTPPPPSSSS((((3333SSSS))))                                                          DDDDLLLLAAAATTTTPPPPSSSS((((3333SSSS))))
  137.  
  138.  
  139.  
  140.      A columnwise scheme is used for solving A*x = b.  The basic algorithm if
  141.      A is lower triangular is
  142.  
  143.           x[1:n] := b[1:n]
  144.           for j = 1, ..., n
  145.                x(j) := x(j) / A(j,j)
  146.                x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
  147.           end
  148.  
  149.      Define bounds on the components of x after j iterations of the loop:
  150.         M(j) = bound on x[1:j]
  151.         G(j) = bound on x[j+1:n]
  152.      Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
  153.  
  154.      Then for iteration j+1 we have
  155.         M(j+1) <= G(j) / | A(j+1,j+1) |
  156.         G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
  157.                <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
  158.  
  159.      where CNORM(j+1) is greater than or equal to the infinity-norm of column
  160.      j+1 of A, not counting the diagonal.  Hence
  161.  
  162.         G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
  163.                      1<=i<=j
  164.      and
  165.  
  166.         |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
  167.                                       1<=i< j
  168.  
  169.      Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTPSV if the
  170.      reciprocal of the largest M(j), j=1,..,n, is larger than
  171.      max(underflow, 1/overflow).
  172.  
  173.      The bound on x(j) is also used to determine when a step in the columnwise
  174.      method can be performed without fear of overflow.  If the computed bound
  175.      is greater than a large constant, x is scaled to prevent overflow, but if
  176.      the bound overflows, x is set to 0, x(j) to 1, and scale to 0, and a
  177.      non-trivial solution to A*x = 0 is found.
  178.  
  179.      Similarly, a row-wise scheme is used to solve A'*x = b.  The basic
  180.      algorithm for A upper triangular is
  181.  
  182.           for j = 1, ..., n
  183.                x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
  184.           end
  185.  
  186.      We simultaneously compute two bounds
  187.           G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
  188.           M(j) = bound on x(i), 1<=i<=j
  189.  
  190.      The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we add
  191.      the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.  Then the
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. DDDDLLLLAAAATTTTPPPPSSSS((((3333SSSS))))                                                          DDDDLLLLAAAATTTTPPPPSSSS((((3333SSSS))))
  203.  
  204.  
  205.  
  206.      bound on x(j) is
  207.  
  208.           M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
  209.  
  210.                <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
  211.                          1<=i<=j
  212.  
  213.      and we can safely call DTPSV if 1/M(n) and 1/G(n) are both greater than
  214.      max(underflow, 1/overflow).
  215.  
  216.  
  217. SSSSEEEEEEEE AAAALLLLSSSSOOOO
  218.      INTRO_LAPACK(3S), INTRO_SCSL(3S)
  219.  
  220.      This man page is available only online.
  221.  
  222.  
  223.  
  224.  
  225.  
  226.  
  227.  
  228.  
  229.  
  230.  
  231.  
  232.  
  233.  
  234.  
  235.  
  236.  
  237.  
  238.  
  239.  
  240.  
  241.  
  242.  
  243.  
  244.  
  245.  
  246.  
  247.  
  248.  
  249.  
  250.  
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  
  261.                                                                         PPPPaaaaggggeeee 4444
  262.  
  263.  
  264.  
  265.